Estudio de lugares de observación óptima para el asteroide 99942 Apophis¶
Librerías¶
In [1]:
"LIBRERÍAS BÁSICAS"
import numpy as np
import pylab as plt
import pandas as pd
import matplotlib.animation as animation
from matplotlib import rc
rc('animation', html='jshtml')
In [2]:
"SPICEYPY"
!pip install spiceypy
import spiceypy as spy
In [3]:
"ASTROPY"
!pip install astroquery
from astroquery.jplhorizons import Horizons
from astropy.time import Time
from astropy.table import Table
from astropy.coordinates import SkyCoord, EarthLocation, AltAz, get_sun
import astropy.units as u
In [4]:
"MAPITAS"
!apt-get install libgeos-3.5.0
!apt-get install libgeos-dev
!pip install https://github.com/matplotlib/basemap/archive/master.zip
!pip install pyproj==1.9.6
from mpl_toolkits.basemap import Basemap
In [5]:
"TQDM"
!pip install tqdm
from tqdm import tqdm
Efemérides¶
In [6]:
"UNIDAD ASTRONÓMICA"
AU = (1*u.au).to(u.m).value
"TIEMPOS"
# Se definen los tiempos desde que está a 12Rt, pasando por su máximo
# acercamiento en 6Rt y volviendo a 12Rt
t_ini = Time('2029-04-13 19:00:00',format='iso')
t_max = Time('2029-04-13 22:05:00',format='iso') #Máximo Acercamiento
# t_end = Time('2029-04-13 23:00:00',format='iso') #Final. Step de 5min
t_end = Time('2029-04-14 06:00:00',format='iso')
timesyy = {'start':t_ini.value, 'stop':t_end.value, 'step':'20m'}
print("(Ignorar advertencia de ERFA si sale, es porque está muy a futuro)")
In [7]:
"EFEMÉRIDES"
#Efemérides de Apofis entre los tiempos estipulados
apophis = Horizons(id='99942',epochs=timesyy)
ephemeris = apophis.ephemerides()
vectors = apophis.vectors()
#Coordenadas Ecuatoriales del Asteroide
ECUcoords = SkyCoord(ephemeris['RA'],ephemeris['DEC']
,ephemeris['delta'],frame='gcrs')
Nt = np.shape(ephemeris)[0]
print(Nt)
In [8]:
# ephemeris.to_pandas()
plt.plot(ephemeris.to_pandas()['datetime_jd']-2.462240e+06,ephemeris.to_pandas()['delta']*23500)
Out[8]:
Coloreando el Mapita¶
In [9]:
def get_h(T,Lon,Lat,RA,Dec,r):
"""
Teniendo las coordenadas ecuatoriales y la ubicación de un observador,
se saca el ángulo de elevación.
"""
G = T.sidereal_time('apparent', 'greenwich').to(u.deg)
LonR = Lon*u.deg + G
#Vector Posición Objeto
Rx = r*np.cos(Dec*u.deg)*np.cos(RA*u.deg)
Ry = r*np.cos(Dec*u.deg)*np.sin(RA*u.deg)
Rz = r*np.sin(Dec*u.deg)
R = np.array([Rx,Ry,Rz])
#Vector Posición Observador
RLx = np.cos(Lat*u.deg)*np.cos(LonR)
RLy = np.cos(Lat*u.deg)*np.sin(LonR)
RLz = np.sin(Lat*u.deg)
RL = np.array([RLx,RLy,RLz])
#Resta de Vectores y Producto Escalar Inverso
R2 = R-RL
H = np.rad2deg(np.arccos(np.dot(RL,R2)/(np.linalg.norm(RL)*np.linalg.norm(R2))))
H = 90 - H
return H
In [10]:
"SACO LA ALTURA DEL ASTEROIDE EN TODAS LAS UBICACIONES (SE DEMORA COMO 2 MINUTOS)"
nlon = 64
nlat = 36
Lons = np.linspace(-180,180,nlon)
Lats = np.linspace(-90,90,nlat)
#Creo el Grid de Latitudes, Longitudes y fechas
Alts = np.zeros((nlon,nlat,Nt))
#Este siempre se demora ratico
for t in tqdm(range(Nt)):
TIME = Time(ephemeris['datetime_jd'][t],format='jd')
RA = ephemeris['RA'][t]
DEC = ephemeris['DEC'][t]
DELTA = ephemeris['delta'][t]*u.au.to(u.earthRad)
for i in range(nlon):
for j in range(nlat):
LON = Lons[i]
LAT = Lats[j]
Alts[i,j,t] = get_h(TIME,LON,LAT,RA,DEC,DELTA)
In [11]:
plt.figure(figsize=(10,6))
plt.contourf(Lons,Lats,Alts[:,:,0].T,levels=np.linspace(-90,90,100),cmap='inferno')
plt.colorbar()
plt.title("Alturas de Apophis en "+str(ephemeris['datetime_str'][0]))
Out[11]:
In [12]:
def filterNeg(arr):
"""
Recibe un array en 3D y le quita los valores negativos
"""
a,b,c = np.shape(arr)
arr1 = np.reshape(arr,a*b*c)
arr1[np.where(arr1 < 0)] = 0
return np.reshape(arr1,(a,b,c))/90
In [13]:
AltsPos = filterNeg(Alts)
plt.figure(figsize=(10,6))
plt.contourf(Lons,Lats,AltsPos[:,:,0].T,levels=np.linspace(0,1,100),cmap='inferno')
plt.colorbar()
plt.title("Alturas Positivas de Apophis en "+str(ephemeris['datetime_str'][0]))
Out[13]:
In [14]:
"ALTURAS DEL SOL (TAMBIÉN SE DEMORA COMO 2 MINUTOS)"
AltsSun = np.zeros((nlon,nlat,Nt))
FilterSun = np.zeros((nlon,nlat,Nt))
#Defino una función de caída para filtrar con la altura del Sol
def sigmoid(x):
return 1 / (1 + np.exp(0.5*(x+10)))
for t in tqdm(range(Nt)):
TIME = Time(ephemeris['datetime_jd'][t],format='jd')
SUN = get_sun(TIME)
Sun_RA = SUN.ra.value
Sun_DEC = SUN.dec.value
Sun_DIST = 1*u.au.to(u.earthRad)
for i in range(nlon):
for j in range(nlat):
LON = Lons[i]
LAT = Lats[j]
Sun_H = get_h(TIME,LON,LAT,
Sun_RA,Sun_DEC,Sun_DIST)
AltsSun[i,j,t] = Sun_H
FilterSun[i,j,t] = sigmoid(Sun_H)
In [15]:
plt.figure(figsize=(10,6))
plt.contourf(Lons,Lats,AltsSun[:,:,0].T,levels=np.linspace(-90,90,100),cmap='inferno')
plt.colorbar()
plt.title("Alturas del Sol en "+str(ephemeris['datetime_str'][0]))
Out[15]:
In [16]:
AltsPosSun = AltsPos * FilterSun
plt.figure(figsize=(10,6))
plt.contourf(Lons,Lats,AltsPosSun[:,:,0].T,levels=np.linspace(0,1,100),cmap='inferno')
plt.colorbar()
plt.title("Alturas de Apophis donde el Sol no Molesta \n en "+str(ephemeris['datetime_str'][0]))
plt.show()
Plot¶
In [17]:
ephemeris.to_pandas().columns
Out[17]:
In [18]:
def get_lon(TIME,RA):
"""
Saca la Longitud que estará bajo el asteroide en cada momento
"""
T = Time(TIME,format='jd')
G = T.sidereal_time('apparent', 'greenwich').to(u.deg)
return RA - G.value
maxlons = get_lon(ephemeris['datetime_jd'],ephemeris['RA'])
maxlats = ephemeris['DEC']
for i in range(len(maxlons)):
if maxlons[i] < -180: maxlons[i] += 360
# for i in range(1,len(maxlons)):
# if maxlons[i] > maxlons[i-1]:
i = -15
maxlons[i] = np.nan
maxlons = np.array(np.insert(maxlons,i+1,[179]))
maxlats = np.array(np.insert(maxlats,i+1,[maxlats[i-1]]))
In [19]:
maxlons[-20:],maxlats[-20:]
Out[19]:
In [20]:
ephemeris['datetime_str'][-12] #Paso Máximo AQUI
Out[20]:
In [21]:
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
t = -1
fig=plt.figure(figsize=(10,8))
ax = plt.axes()
ax.clear()
plot = ax.contourf(Lons,Lats,AltsPosSun[:,:,t].T*90,levels=np.linspace(-0.1,1,100)*90,cmap='Blues')
m = Basemap(ax=ax,projection='cyl', lon_0 = 0, lat_0 = 0)
m.drawcoastlines()
m.drawparallels(np.arange(-90,90,30),labels=[1,1,0,1])
m.drawmeridians(np.arange(-180,180,30),labels=[1,1,0,1])
ax.set_title("Trayectoria y Elevaciones para el Asteroide 99942 Apophis \n"+str(ephemeris['datetime_str'][t]))
ax.plot(maxlons,maxlats,'w-',lw=5)
ax.plot(maxlons[t],maxlats[t],'wo',ms=12)
ax.plot(maxlons,maxlats,'b-',lw=2)
ax.plot(maxlons[t],maxlats[t],'bo',ms=9)
axins = inset_axes(ax,
width="5%", # width = 5% of parent_bbox width
height="100%", # height : 50%
loc='lower left',
bbox_to_anchor=(1.07, 0., 1, 1),
bbox_transform=ax.transAxes,
borderpad=0,
)
cbar = fig.colorbar(plot, cax=axins)
cbar.ax.get_yaxis().labelpad = 15
cbar.ax.set_ylabel('Elevación sobre el Horizonte [°]', rotation=270)
plt.show()
In [22]:
def frame(t):
ax.clear()
plot = ax.contourf(Lons,Lats,AltsPosSun[:,:,t].T*90,levels=np.linspace(-0.1,1,100)*90,cmap='Blues')
m = Basemap(ax=ax,projection='cyl', lon_0 = 0, lat_0 = 0)
m.drawcoastlines()
m.drawparallels(np.arange(-90,90,30),labels=[1,1,0,1])
m.drawmeridians(np.arange(-180,180,30),labels=[1,1,0,1])
ax.set_title("Mapa de Contorno de Elevaciones y Trayectoria para el Asteroide 99942 Apophis \n"+str(ephemeris['datetime_str'][t]))
ax.plot(maxlons,maxlats,'w-',lw=5)
ax.plot(maxlons[t],maxlats[t],'wo',ms=12)
ax.plot(maxlons,maxlats,'b-',lw=2)
ax.plot(maxlons[t],maxlats[t],'bo',ms=9)
axins = inset_axes(ax,
width="5%", # width = 5% of parent_bbox width
height="100%", # height : 50%
loc='lower left',
bbox_to_anchor=(1.07, 0., 1, 1),
bbox_transform=ax.transAxes,
borderpad=0,
)
cbar = fig.colorbar(plot, cax=axins)
cbar.ax.get_yaxis().labelpad = 15
cbar.ax.set_ylabel('Elevación sobre el Horizonte [°]', rotation=270)
return plot
Animación¶
In [23]:
fig=plt.figure(figsize=(10,8))
ax = plt.axes()
anim = animation.FuncAnimation(fig, frame, frames=Nt, blit=False, repeat=True)
In [24]:
Writer = animation.writers['ffmpeg']
writer = Writer(fps=5, metadata=dict(artist='Agustín-Vallejo'), bitrate=1800)
anim.save("Apophis3.mp4",writer=writer)
from google.colab import files
files.download("Apophis3.mp4")
In [25]:
anim
Out[25]:
In [26]:
from google.colab import drive
drive.mount('/content/drive')
FAST¶
In [ ]:
"UNIDAD ASTRONÓMICA"
AU = (1*u.au).to(u.m).value
"TIEMPOS"
# Se definen los tiempos desde que está a 12Rt, pasando por su máximo
# acercamiento en 6Rt y volviendo a 12Rt
t_ini = Time('2029-04-13 19:00:00',format='iso')
t_max = Time('2029-04-13 22:05:00',format='iso') #Máximo Acercamiento
t_end = Time('2029-04-14 12:00:00',format='iso')
timesyy = {'start':t_ini.value, 'stop':t_end.value, 'step':'15m'}
print("(Ignorar advertencia de ERFA si sale, es porque está muy a futuro)")
In [ ]:
"EFEMÉRIDES"
#Efemérides de Apofis entre los tiempos estipulados
arecibo = { 'lon': 293.24692,
'lat': 0.949577,
'elevation': +0.312734}
fast = { 'lon': 106.85667,
'lat': 25.653055,
'elevation': +0.312734}
apophis = Horizons(id='99942',location=fast,epochs=timesyy)
ephemeris = apophis.ephemerides()
apophis.vectors()
In [ ]:
ephemeris.to_pandas().columns
Out[ ]:
In [ ]:
ephemeris['horas'] = 24*(ephemeris['datetime_jd']-ephemeris['datetime_jd'][0])
In [ ]:
H1 = ephemeris['horas'][np.where(ephemeris['EL'] >= 50)][0]
H2 = ephemeris['horas'][np.where(ephemeris['EL'] >= 50)][-1]
print(H1,H2)
In [ ]:
ephemeris[['datetime_str','RA','DEC','EL']][np.where(ephemeris['EL'] >= 49)].to_pandas()
Out[ ]:
In [ ]:
fig,ax = plt.subplots(figsize=(8,5))
ax.set_title("Efemérides de 99942 Apophis desde FAST (25.6°N - 106.8°E)")
ax.plot(ephemeris['horas'],ephemeris['EL'],label='Elevación [°]',lw=3,color='b')
# plt.plot(ephemeris.to_pandas()['horas'],ephemeris.to_pandas()['DEC'],label='Declinación [°]',lw=3)
# plt.plot(ephemeris.to_pandas()['horas'],ephemeris.to_pandas()['elong'],label='Elongación [°]',lw=3)
ax.fill_between([H1-0.2,H2],0,100,alpha=0.2,color='b')
ax.text(8.5,45,"Apophis Visible \n desde FAST",color='b')
# ax.hlines(50,-10,50,'b',linestyles='--')
ax.set_ylabel("Elevación del Apophis (°)",color='b',size=14)
plt.xlabel("UT el 2029-04-14",size=15)
plt.xticks([0,5,10,15],['-5:00','0:00','5:00','10:00'],size=15)
plt.yticks(size=12)
plt.ylim([0,90])
plt.xlim([0,ephemeris['horas'][-1]])
ax.grid()
ax2=ax.twinx()
ax2.plot(ephemeris['horas'],ephemeris['delta']*23500,label='Distancia [Rt]',lw=3,color='r')
ax.plot(ephemeris['horas'],ephemeris['delta']*0-1,label='Distancia [Rt]',lw=3,color='r')
ax2.set_ylabel("Distancia al Apophis ($R_t$)",color='r',size=14)
ax.legend(loc=2)
plt.show()
